Doubly degenerate orbital system in honeycomb lattice: 
implication of orbital state in layered iron oxide 
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We study a doubly-degenerate orbital model on a honeycomb lattice. This is a model for orbital states in 
multiferroic layered iron oxides. The classical and quantum models are analyzed by spin-wave approximation, 
Monte-Carlo simulation and Lanczos method. A macroscopic number of degeneracy exists in the classical 
ground state. In the classical model, a peak in the specific heat appears at a temperature which is much lower 
than the mean-field ordering one. Below this temperature, angle of orbital pseudo-spin is fixed, but conventional 
orbital orders are not suggested. The degeneracy in the ground state is partially lifted by thermal fluctuation. 
We suggest a role of zero-dimensional fluctuation in hexagons on a low-temperature orbital structure. Lifting 
of the degeneracy also occurs at zero temperature due to the quantum zero-point fluctuation. We show that the 
ground-state wave function is well represented by a linear combination of the states where a honeycomb lattice 
is covered by nearest-neighboring pairs of orbitals with the minimum bond energy. 

PACS numbers: 75.30.-m, 71.10.-w, 75.10.Jm 






I 

O 

o 



> 

en 

00 

o 
^' 

O 
00 

o 



X 



I. INTRODCUTION 

Orbital degree of freedom and its interplay with spin and 
charge degrees are one of the recent attractive themes in con- 
densed matter physics.-i^ Orbital represents an anisotropic 
shape of the electronic wave function. In a molecule, this de- 
gree of freedom is quenched by the Jahn-Teller effect, and/or 
a formation of the chemical bond along a specific bond di- 
rection. On the contrary, in a solid crystal, some equivalent 
bonds coexist. One alignment of orbitals dose not fully satisfy 
the minimum-energy configuration for all equivalent bonds. 
This is a certain kind of frustration subsisting intrinsically in 
a solid crystal with orbital degeneracy. This frustrating and 
directional character for the orbital provides a wide variety of 
exotic phenomena in transition-metal compounds near a Mott 
insulating state. 

For orbital degenerate systems under strong electron corre- 
lation, a number of theoretical investigations have been done 
for more than one decade. One of the well known and ex- 
amined orbital models is the so-called three-dimensional eg 
orbital modeli^i^ This is proposed as a model for orbital state 
in LaMnOs and KCUF3 with the perovskite crystal structure. 
The doubly degenerate Cg orbitals, d^2_^,2 and dj^-i_p., are rep- 
resented by the pseudo-spin (PS) operator T with magnitude 
of 1/2 and are located on a simple cubic lattice. The model 
Hamiltonian is given by 



^«.=-^E <^^e,. 



^'^, 



(1) 



Here, a vector Ctj for 77 = [x, y, z) connects the nearest neigh- 
boring (NN) sites, and x^ is a linear combination of the PS 
operator defined by t/' = — sin(2;r«,] /3 ) T^ + cos {iKn-q /3 ) T^ 
with a factor (nv,nv,«z) = (li2,3). This model is derived 
by the perturbational procedure from the extended Hubbard 
Hamiltonian with neglecting spin degree of freedom. The r\ 
dependence of the interaction implies the frustrating and di- 
rectional character. As seen in frustrated magnets, there is 



a macroscopic number of degeneracy in the classical ground 
state. This degeneracy is lifted by thermal fluctuation in finite 
temperatures and by quantum zero-point fluctuationj ^-^i^-^ As 
results, a staggered-type long-range orbital order is realized. 

Doubly-degenerate orbital model on a honeycomb lattice, 
studied in the present paper, is one of the orbital models with 
the frustrating and directional interaction. Orbital degree of 
freedom represented by the PS operator is located on a two- 
dimensional honeycomb lattice [see Fig.[T]. An explicit form 
of the Hamiltonian is given in Eq. (|4]i, which is introduced in 
more detail in Sect. Ill Al This model looks similar to the eg 
orbital model in Eq. ([U; the NN three-bond directions in a 
honeycomb lattice, a, /3, and 7, correspond to the Cartesian 
coordinates in a cubic lattice. Thus, a similar kind of frustrat- 
ing character for orbital configuration is expected. However, 
in general, stability of an orbital state is extremely sensitive to 
symmetry and dimension of a crystal lattice. It is nontrivial 
whether a conventional long-range order is realized or not, in 
the same type of interaction, but in the different crystal lattice. 
From a viewpoint of substantial materials, the honeycomb- 
lattice orbital model is proposed as an orbital model in a mul- 
tiferroic layered iron oxide R¥e20/\ (R=L\x, Y, Yb) i''-'^ This 
is a mixed valence compound where equal amount of Fe^+ 
and Fe^+ coexists in a pair of triangular lattice ji^ii^*i^iiii^ A 
Fe^+ ion with cfi configuration has the doubly degenerate or- 
bital degree of freedom. In the low-temperature charge and 
spin ordered phase, a Fe^+ sublattice forms a honeycomb lat- 
tice, and the orbital state is mapped onto a honeycomb lattice 
model. This will be introduced later in more detail. From dif- 
ferent view point, this orbital model is proposed recently in 
study of the optical lattice i'V^-^' 

In this paper, we study the ground-state and finite- 
temperature properties in the doubly-degenerate orbital model 
on a honeycomb lattice. We analyze the classical and quan- 
tum models by the Monte-Carlo (MC) and Lanczos methods, 
respectively, as well as the spin-wave approximation. There 
is a number of the degenerate classical ground states as well 
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FIG. 1: (a) A honeycomb lattice structure and sublattices A and B. 
Bold arrows represent vectors connecting NN two sites, (b) Brillouin 
zone and reciprocal lattice vectors for a honeycomb lattice. 
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FIG. 2: Pseudo-spin operator T, and its projection components T^ 
along the three-bond directions. 



as the eg orbital model. In the classical model, at a certain 
temperature which is much lower than the mean-field order- 
ing temperature, a peak in the specific heat appears. Below 
this temperature, the PS angles are fixed at one of 7Zn/3 with 
an integer number «. The degeneracy is partially lifted below 
this temperature due to thermal fluctuation, but the conven- 
tional long-range orders are not suggested from the two-body 
correlation functions for PS. This degeneracy is also lifted by 
the quantum zero-point fluctuation. The ground-state wave- 
function is well reproduced by a linear combination of the 
states that a honeycomb lattice is covered by dimer pairs of 
the NN PS configurations which satisfy the minimum bond 
energy. 

In Sect. Ill A! we define the Hamiltonian of the honeycomb 
lattice orbital model, and introduce implication of the orbital 
state in layered iron oxides. Results in the classical and quan- 
tum models are presented in Sects. |III] and |IV| respectively. 
Section rvlis devoted to the discussion and summary. Prelimi- 
nary results have been published in Refs. |Tl] and[i2t Relation 
to the layered iron oxides is briefly introduced in Ref.[l3i 



II. MODEL 

A. Model Hamiltonian 

We start with the model Hamiltonian for the doubly degen- 
erate orbitals, denoted by a and b, defined in a honeycomb 
lattice. This is represented by the pseudo-spin operator with 
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FIG. 3 : Eigen values of the orbital interaction 7(k) in the momentum 
space. 



magnitude of 1/2: 



T; = :^L4.-^«'4'. 



(2) 



where djts is the electron annihilation operator with orbital f(= 
a,b), spin ^(=1,4) at site /, and a are the Pauli matrices. For 
the three-kinds of NN bonds, Tj = (a,j8,7), in a honeycomb 
lattice (see Fig.[rii, we introduce new PS operator as 



iTttlr 



T; 



litrir 



T:\ 



(3) 



A numerical factor «,j is defined as («Q;,«jg,«y) = (0,1,2). 
When we define the pseudo-spin coordinate as shown in 
Fig. |2] the operator T^ represents a projection component of 
T, on the Tj bond direction. The model Hamiltonian studied 
in the present paper is. 



J^- 



ieA 



''i ''i+ea 



Tf V 



S- '■i+i 



(4) 



where e^ is a vector connecting the NN sites along the direc- 
tion T], i^jgA represents a sum of sites on the sublattiece A 
[see Fig. [Ua)], and J is the exchange constant. Although J 
is defined to be positive, its sign is gauged away by rotating 
PS's on the A sublattice with respect to T' . This Hamiltonian 
is rewritten as a following simple form 



Jf: 



^ ieA.,ri 



i'-tI 



- ieA 



rx2 



rz2\ 



(5) 



The second term is — 3 JN/ 16, when Tj is a two-dimensional 
classical spin, and is —3JN/di in the quantum-spin case. A 
total number of sites is A^. This model is proposed as a orbital 
state for the layered iron oxide , ^''^^ as explained in Sect. IIIBI 
in more detail, and is also recently proposed in study of the 
optical lattice i'VQ-^' A similar orbital model in a honeycomb 
lattice termed the Kitaev model is recently well examined i^^'^^ 
Here three components of the PS operator, 7j' with / = {x,y,z), 
instead of t/ in this model, are concerned in the interactions 
along the a, j8 and y directions. 



Before going to detailed analyses of the Hamiltonian, we 
briefly introduce a character in this model. Let introduce the 
Fourier transformation for the orbital PS, 



Tc(k) 



1 



ET,e''^"', 



(6) 



for the sublattice C(=A, B). The Hamiltonian Eq. (HJi is repre- 
sented in the momentum space )'^''^ shown in Fig.[TJb), as 



=^=i/'(-k)/(k)v/(k) 



(7) 



We introduce a four-component vector defined as 

V/(k) = K(k),r^(k),7^'(k),7B^(k)] , (8) 

and a 4x4 matrix -/(k). We obtain 

the eigen values of y(k) are ±37/4 and 
±y[3+2cosk-a + 2cosk-b + 2cosk- (a-b)]'''^/4 where 
a and b are the primitive translation vectors defined in Fig. [T] 
Numerical plot of 7(k) is presented in Fig. |3] The lowest 
eigen value is a momentum independent flat band of —37/4. 
That is, the effective dimensionality for the lowest state is 
zero, and, in the classical ground state, stable orbital struc- 
tures are not determined uniquely due to large fluctuation. 
The second eigen value touches the lowest band at the point 
F. 

Compare the present model with the e^, orbital model in a 
simple cubic lattice. The Cg orbital model defined in Eq. ([TJ 
shows a similar form with the present honeycomb lattice 
model in Eq. (|4]i, when a, j3 and 7 are replaced by the Carte- 
sian coordinates x, y and z- The momentum representation 
of the orbital interaction is given by J{k) ~ ±J[3 + cosk^a + 
cos kya + cos k^aY'^ where {kx,ky,k.) are defined in the Bril- 
louin zone for a simple cubic latticei^ Dispersion relation of 
/(k) is flat along (tt, tt, k) — (0, tt, k) and other equivalent di- 
rections. Due to the flat dispersions, there is a macroscopic 
number of degeneracy in the classical ground state. However, 
this degeneracy is lifted by thermal and quantum fluctuations, 
and a staggered long-range orbital order is realizedj^ This 
is the so-called order-by-fluctuation mechanism. The long- 
range order in the classical model is confirmed by the Monte- 
Carlo simulation; the two-body correlation function for PS at 
momentum k = (tt, tt, k) starts to increase around T = 0.177, 
and is saturated at its maximum value in the low temperature 
limit [see inset in Fig. fTTT c'll.i^^ 

B. Implication of layered iron oxide 

In this section, we introduce the honeycomb lattice orbital 
model defined in Eq. ^ as an orbital model for multiferroic 
layered-iron oxides RFcxOa- This is known as a multiferroic 
material driven by electronic charge and spin degrees of free- 
dom. Electric and magnetic properties in /?Fe204 are domi- 
nated by Fe 3d electrons in a pair of triangular-lattice planes 
stacked along the c axis, which is termed the W-layer [see 
Fig. |4ja)]. A Fe ion in the W-layer is five fold coordinate 
with a local symmetry of Dm- The five 3d orbitals under the 




FIG. 4: (a) A pair of triangular planes termed the W-layer, and (b) 
three Fe-O bond directions in a triangular lattice in WFe204. 



crystalline field split into two sets of the doubly degenerate or- 
bitals, {dxy,d^2_^,2} with the symmetry E', and {dy-, d-^} with 
E", and the di,^2_,.2 orbital with A'. We obtained by the crys- 
talline field calculation that the E' orbital is the lowest. Since a 
nominal valence of the Fe ions is 2.5h-, equal amount of Fe^+ 
(d^) and Fe^+ (d^) coexists. The five 3d orbitals are singly 
occupied in Fe^+, and one of the degenerate lowest orbitals in 
Fe^+ is doubly occupied. Thus, Fe^+ has the doubly degen- 
erate orbital degree of freedom. This is represented by the PS 
operator defined in Eq. (|2| where t takes dxy and d^2_^,2- It is 
convenient to introduce the three two-dimensional coordinates 
(77(, 77,') with 77' = (a',)3', /) where the 77^ axis is parallel to 
one of the NN Fe-O bonds as shown in Fig.Ub). We define, 
in these coordinates, linear combinations of the orbital opera- 
tors: 



d. lo '1 
ir\/-ri,.-s 



r]'X.s 



d, 



Au • A-Tt 

— sm^njj/, cos^njj/ 



dix^-,2, 



^ix\s 



(9) 



with a numerical factor {n(j^i^nai,ny) =- (0,1,2). In the NN 
Fe-O bond along the r\' axis, the d 1212 and O 2p orbitals 
form the a bond. We redefine the PS operators. 



n' / 2?r \ _ ( 1% . 

T '- = COS I — n,,/ 1 7;.- + sm ( —n^i \ T; . 



(10) 



One hole occupied state in the d'l'i (d 

'Ix ^y 



v'.ni' 



orbital at site / 



is the eigen state of T^ . 

Interaction between the orbitals is constructed from the 
electronic processes in a W-layer. The model Hamiltonian 
in low energy spin, charge and orbital states is derived from 
the extended pd model by the perturbational procedure. The 
obtained Hamiltonian consists of the long-range Coulomb in- 
teractions between charges and the exchange interactions be- 
tween NN spins and orbitals. We analyze numerically the 
Hamiltonian by the classical MC method. Details were pre- 
sented in Refs. [ij and [Tj. Obtained charge and spin or- 
dered structure is shown in Fig. |5] which is consistent with 
the electron and neutron diffraction experiments , 1 5, 1 6, 1 7, 1 8 ^ 
charge imbalance of Fe^+ and Fe^+ is realized between the 
triangular-lattice planes. That is, the electric dipole moment is 







FIG. 6: Pseudospin configurations in the ground state. Arrows rep- 
resent directions of PS's and bold bars are for the projection compo- 

FIG. 5 : Schematic picture of the charge and spin structures in 2Fe2+ - "^nts T^ along the bond direction. A PS configuration, obtained by a 

Fe3+ plane (right) and in Fe2+-2Fe3+ plane (left) for ^Fe204. Filled uniform rotation of all PS's in right figure, is also in the ground state. 

and open circles represent Fe^+ and Fe^+, respectively. At sites 

surrounded by dotted circles, spin directions are not uniquely de- 

termined due to frustration. A. Orbital structure at ground state 



caused by the charge order without inversion symmetry J^ In 
the spin structure shown in Fig.|5j Fe^+ ions, which have the 
orbital degree of freedom, are surrounded by NN Fe^+ in the 
Fe^+ -2Fe^+ plane, and these form a honeycomb lattice in the 
2Fe^+ -Fe^+ plane. The superexchange interactions in three 
Fe^+ -Fe^+ bonds connecting Fe^+ at site i is proportional to 

Y^rii T^ . This is because the orbital is only active in Fe^+ , 
and spin configurations in the three bonds are equivalent. It 
is easily shown from Eq. dTol ) that this is zero. Therefore, the 
orbital degree of freedom in the charge and spin ordered phase 
is described by the Hamiltonian in a honeycomb lattice in the 
2Fe^+-Fe^+ plane. 



J^' 



-J'L 

i'GA 



Tf T 



i a! 



(11) 



Orbital structure in the classical ground state is obtained 
from the Hamiltonian in Eq. (|5]l. The ground state energy is 
—37/16, when the PS's satisfy the following condition in all 

NN bonds; '^1^°!^' 

^;'-4e,- (13) 

This relation implies that the projection components of PS's 
are equal with each other for all NN bonds. There is a macro- 
scopic number of orbital structures which satisfy this condi- 
tion. Two of them are shown in Fig. |6l In particular, uniform 
orbital alignments with any PS angles are in the ground state 
configurations. This kind of rotational symmetry is not ex- 
pected from the Hamiltonian where any continuous symme- 
tries do not exist in the PS space. But this is consistent with 
the momentum representation of the orbital interaction, 7(k); 
the second-lowest band in y(k) touches the lowest one at the 
point r as shown in Fig. [3] 



The exchange constant /(> 0) is given by the intra-site 
Coulomb interactions and the hopping integrals. Then, we 
introduce the unitary transformation. 



'^-A-'W-^w 



(12) 



which rotates PS's on sublattice A(B) by angle 7r/6 (5;r/6) 
with respect to the T^' axis. We show that U^^J^/f'U is identi- 
cal to J^^ in Eq. (|4]i where J corresponds to /. In addition to 
the exchange interaction described by this Hamiltonian, there 
may be some other factors which couple with orbital degree of 
freedom. However, this Hamiltonian is expected to provide a 
starting point to examine the low temperature orbital structure 
in layered iron oxides. 



III. CLASSICAL ORBITAL STATE 

In this section, we treat the orbital pseudo-spin T, as a clas- 
sical two-dimensional vector with an amplitude of 1/2. 



B. Spin wave analyses 

At the first step, among the degenerate uniform configura- 
tions at zero temperature, we turn up stable states in finite tem- 
peratures by using the spin-wave approximation!^ We define 
the PS angle as 0,- = — tan^'(7jY7j'^), and denote an angle in 
the uniform configuration by 0*. A deviation from 0* at site 
/ is represented by ^,(= 0,-0*). Within the second order of 
^i, the spin-wave Hamiltonian is obtained as 



/ 



^SW = T L ?r,(0*)(C,-C,4 



(14) 



■ iGA 77 



where ^,,(0*) = \sm^[9* + (27r«,,)/3]. 


By in- 


troducing the Fourier transform of Q 


defined by 



Ck == (A^/2)"'/^I,ec «*■'■' C; for sublattice C(=A,B), the 
Hamiltonian is rewritten in a momentum space as 



J 



=^sw = :ri:?r,(r)5:|Ck^-Ck^. 



B „—ik-en 



(15) 



Then, we calculate the partition function for the PS fluctuation 
around 0*. By introducing the two-dimensional polar coordi- 
nates defined by 1^^ — |Ck'|e"''k for C=A and B, the partition 
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FIG. 7: Contour map of the function /(0*,k) in the Brillouin zone 
for e* = in (a), and that for 0* = :;r/6 in (b). 
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FIG. 8: Left: PS configuration for 0* = 0. Right: configuration 
obtained by ±59 rotations of PS's in each zigzag chain. 



function is obtained as 



X exp 



I? 



m-\^^w^'^e 



.BU'A-Pk^-'ke^ 



(:<5) 



where A(> 0) is the Jacobian, j5 is the inverse temperature, 
^(Pk = "Pk ~ "Pk^' ^nd njj represents a product of k in a half of 
the first Brillouin zone. At low temperature, the upper limits 
in the integrals for | ^j^ | and | ^^ | are safely taken to be infinity. 
By integrating out a variable |Ci^|, we obtain the following 
expression for the free energy, 

F{e*) = -liogA-^iog 

- il/(0*.k), 

P k 



iPJ)' 



(17) 



with 



/(0*,k) = log 



d^ I d(p 



[Lnqvmil-l^e'fe 



-!k-e„ 121 



(18) 



where ^j^ represents a sum of k in a half of the first Brillouin 
zone. 

We numerically calculate /(0*,k). Contour maps of 
/(0*,k)for0* = and ;r/6 are presented in Fig. |7] Results in 



<fc. 




FIG. 9: Apart of the free energy F(0*) as a function of the PS angle 
0* obtained in the spin wave approximation. 



other 6* ~ 2nn/6 and (2« + 1 ) n/6 with integer n are obtained 
by considering the Cg symmetry in /(0*,k). This symmetry 
is attributed to the fact that the Hamiltonian is invariant under 
(i) the inversion with respect to PS, and (ii) a combined oper- 
ation of the C3 rotation for PS and that for the crystal lattice. 
In f{9* = 0,k), a divergent behavior appears along the Gq 
(horizontal) axis. This originates from a number of low-lying 
PS configurations from the 0* = state, explained as follows. 
Start with the PS configuration with 0* = shown in Fig. [S] 
and focus on zigzag chains running along the b (vertical) axis. 
Rotate PS's by angle +50 or —50, where |50| is taken to be 
uniform and their signs are chosen independently for the each 
zigzag chain. One example is shown in Fig. [8] This rotation 
does not change the energy, since the condition in Eq. ( fT3b is 
still satisfied in all NN bonds. On the contrary, in 0* = 7r/6, 
a divergent behavior in /(0*,k) is only seen at the point F. 
This corresponds to a uniform PS rotation. By integrating out 
the momentum k for /(0*,k), we obtain the 0* dependence 
of the free energy. We present, in Fig. |9l a part of the free 
energy defined by 



F{e*) = -^Xf^e*,k). 



N' 



(19) 



Because of the one-dimensional fluctuation in /( 0* , k), F ( 0* ) 
takes its minima at six angles of 0* = n7z/3. An analytical 
form is given as F(0* =«7r/3) = - log( 16/3) — (1/2) log ttw 
—2.246. Among the continuous uniform states, these six 
states are stabilized selectively by thermal fluctuation. 



C. Monte Carlo simulation 

In the previous section, we assume the uniform PS con- 
figurations and show lifting of the continuous degeneracy by 
thermal fluctuation within the spin-wave scheme. Now we 
take off this restriction and show the results obtained by the 
MC simulation. Because of a limited system size in the MC 
calculation, both the spin wave and MC methods provide us 
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FIG. 10: Specific heat calculated in several clu.ster sizes. Low tem- 
perature data are enlarged in (b). Inset in (b) shows a peak position 
To in the specific heat as function of 1/A'. 
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FIG. 11: Correlation functions 5'"'(k) for several momentum q. 
Cluster sizes are 2 x 2 x 2 in (a), 2 x 3 x 3 in (b), and 2x6x6 
in (c). Inset in (c) shows correlation function 5(k) = 4N^ Y,ij{Ti ' 
Tj) £'"'■(■"'"''>' at k = {n,n,n) calculated in the eg orbital model. A 
cubic cluster with 1 S"' sites is adoptedi2ii2 



complemented information with each other To avoid a trap of 
the simulation in local minima, we adopt the multi-canonical 
MC technique. The energy distribution functions are obtained 
by the histogram method-- and the CFP ona^. In most of 
the simulation, 1x10^ MC steps are used to produce the en- 
ergy histogram, and 2x 10^ MC steps are for the calculation. 
Statistical averages and errors are obtained by 20 times sim- 
ulations. Except for the results in Fig. [TSlb). error bars are 
enough small and are not plotted in the figures. We adopt a 
cluster of 2 X L X L(= A^) sites with L = 2 -- 24. 

First we present, in Fig. [10] temperature dependence of 
the specific heat C{T) for several system sizes. As seen in 
Fig. [Tol a), over all behavior does not show size dependence. 
There is a shoulder around 0.17 and a sharp peak around 
0.0057 — 0.017 which depends on system size. Result in a 
2x5x7 size cluster is almost identical with that in 2 x 6 x 6; 
a shape of the cluster is not essential. A magnification of C(r) 
in a low temperature region is presented in Fig. [TOlb). With 
increasing a system size, the peak shifts to a lower tempera- 
ture side and becomes sharp. The peak position is denoted as 
To from now on. As shown in the inset of Fig. [TOl b). To ap- 
proaches a finite value about 0.00647 in the thermodynamic 
limit. It is worth noting that this value of To is much smaller 
than the mean-field ordering temperature 37/8. At zero tem- 
perature limit, C{T) takes about 0.5 corresponding to one de- 



gree of freedom per site, i.e. the two-dimensional PS angle. 

To elucidate the PS structure below 7b, we calculate the 
correlation functions for PS defined by 



5''"(k) 



4 



YjjM 



,'k-(r,- 



(20) 



where / and m take x and z, and r, is a position of site /. The 
maximum value of the functions is one. The z component 
of the correlation functions ^"(k) for several system sizes are 
presented in Fig.[TT] We calculate 5'^^(k)'s for all possible mo- 
menta k in a cluster In a 2 x 2 x 2 cluster, 5^^(k) at k = (0, 0) 
takes about 0.3 in low temperatures. However, with increas- 
ing A^, the values of 5''"(k) decrease rapidly, and in a 2 x 6 x 6 
cluster, all 5''"'(k)'s are less than 3% of their maximum value. 
Other components, ^"(k) and 5''''^(k), are similar to ^^(k). 
We conclude that, below 7b, there are no conventional long- 
range order corresponding to the correlation functions given 
in Eq. ( [20l i. This is not trivial for the present model where 
the Mermin- Wagner's theorem is not applicable. The present 
results are in contrast to those in the eg orbital model; the PS 
correlation function at k = (tt, tt, n) starts to increase around 
0.177, and approaches its maximum value at the low temper- 
ature limit, as shown in the inset of Fig. [TTT c). 

Here we propose a physical parameter q for the PS angle 0, 
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FIG. 12: (a) Correlation function of a variable q for the PS angle. 
Insets show the q=\ and — 1 PS configurations, (b) Scaling analyses 
for the correlation length of qi = cos 30/. 



defined by 



^=-^COs30;. 



A^ 



(21) 



Because the Hamiltonian is invariant under the inversion of 
all Ps's, (q) is zero in a disordered phase. When the angle 0, 
takes oneof the three angles 2«7r/3 [{2n + l)n/3],q= 1 (—1) 
[see inset of Fig. [T2l" a)l. In Fig. [121 we plot the temperature 
dependence of the correlation function of q defined by 

(22) 

This starts to increase around To and is saturated to the max- 
imum value at the low temperature limit. With increasing the 
system size A^, Q abruptly increases around Tq. We show, in 
Fig. fT2l b). the finite-size scahng plot for the correlation length 
of qi = cos 3 0, defined by^ 



1 dG{k) 



G(0) £/|k|2 



|k|=0 



where we define the correlation function of qi as 

1 



G{k)^-^j:{q.qj)e''^-^^--'^^K 



N- 



(23) 



(24) 



As shown in this figure, £,/L in several system sizes are 
scaled by the scaling function {T — 7b)L''*'/7b within er- 
ror bars. Here we obtain the transition temperature Tq = 



-0.5 " 




FIG. 13: Temperature dependence of congelation functions of qi and 
specific heat. Broken, dotted and dash-dotted lines are for the corre- 
lation functions between the NN, 2nd NN and 3rd NN sites, respec- 
tively. 



0.0067 ± 0.0007 and v = 0.72 ±0,04. These results imply 
that, at low temperature below Tq, the PS angle at each site 
takes one of the three angles 2mt/3 (cos30 = 1), or one of 
(2n+ l);r/3 (cos 30 = -1). When the qi = 1 and -1 states 
are randomly distributed in a lattice, the correlation function 
Q should be zero. It is found from the snapshot of the MC 
simulation that the three qi =1 states, or the three —1 states, 
coexist below Tq. From the view point of the PS angle, a 
shoulder structure in C(r) around T// = 0.1 shown in Fig.fTOl 
coiTesponds to development of the short range coiTelation. In 
Fig. [13] we show the short-range correlation functions of qi 
defined by 



Q('n) 



1 



j('n)A? 



Y.'^lilj) 



(25) 



{'■;■) 



where G^'"' with m =1,2 and 3 are the correlations between 
NN, the next NN and the 3rd NN sites, respectively. A numer- 
ical factor z^™' is a number of the neighboring pairs, and ^', . ■■, 
represents a sum of the pairs. It is clearly shown in Fig. [TS] 
that a shoulder of C(r) corresponds to development of G^''. 

Stability of the q = ±.\ states is attributed to the low-lying 
excited states around the q = ±.\ states. Consider one of the 
q = \ states shown in Fig. [Ha), and local PS fluctuations 
in a certain NN bond in this state. There are two ways for 
fluctuation where the condition in Eq. ( [T3] ) is satisfied in this 
bond. That is, the two kinds of excited states appear thermally 
with the same probability. Situation is different away from 
the q= \ configuration. Consider one of the q ^ ±\ states 
shown in Fig. [T4|b). There is only one way for fluctuation 
where the condition in Eq. ( [T3] l is satisfied. This high density 
of the low-lying fluctuations around q = ±.\ states contributes 
to the entropy gain and stabilizes the q = ±1 states at finite 
temperature. ■- 

As shown above, we have found that, below Tq, the PS an- 
gle at each site is fixed at one of the three angles 2n7r/3 or 
one of the three [In + 1 );r/3. Within the present calculations, 
we do not insist whether all ^ = ±1 states are realized equiv- 
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FIG. 14: Schematic view of tiie q - 
low-lying fluctuations. 
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FIG. 15: Energy correction AE(6*) from the zero-point vibration as 
a function of the orbital angle 0* . 



alently or not. We will discuss this point in Sect. [V]in more 
detail with supplementary calculations. 



IV. QUANTUM ORBITAL STATE 

In this section, we analyze the Hamiltonian in Eq. (|4]i where 
the PS is treated as a quantum spin operator with a magnitude 
of5=l/2. 



A. Spin wave analyses 

To elucidate roles of the quantum fluctuation on the stable 
orbital state at zero temperature, we start from, for simplicity, 
the uniform orbital state with an PS angle 0*. The Hamilto- 
nian is analyzed by using the Holstein-Primakoff transforma- 
tion4 We utilize the rotating frame given by the unitary trans- 
formation with respect to the T^ axis, and introduce the two 
Holstein-Primakoff bosons, au and b^, for the two sublattices 



in the honeycomb lattice. The Hamiltonian up to the second 
order of the boson operator is given as 



^^s 



sw 



Is'-jN+lsjy 



ala^ + blpk 



t;,t 



2^ii(0*)«k&-k+7-k(0*K&:k 

+ ll(0*)ak^i + 7-k(0*K^k 



where 7k(0*) is the structure factor defined by 



^(0*) = |^sin2fr-^„, 



..-ik-en 



(26) 



(27) 



with a numerical factor {iia,n^,ny) ~ (0, 1,2). The first term 
in Eq. ( |26] ) is the zero-th order energy, denoted by £0. which 
corresponds to the second term of Eq. (|5]) in the classical spin 
case. This energy is independent of the angle 0*, as men- 
tioned previously. By applying the Bogoliubov transforma- 
tion, we diagonalize the second term and obtain. 



jot m 



^w=£o+A£(0*)+ ^ <(0*)4"v; 



k ' 



(28) 



k/=± 



„(±) 



where we introduce the boson (orbiton) operators cj^ and 
their energy dispersions 



ai^\e*) = \sj^\±\UQ* 



(29) 



In the case of Q* ~ nK/3, these dispersions show an 
one-dimensional character; for example, O^ [9* = 0) = 
(357/2) -y/l ±cos(^yfl/2) which is independent of kx, where 
we define fc, = k • Go and fcj. = k • (Ga + 2G/,)/\/3. The sec- 
ond term in Eq. ( |28] ) corresponds to the correction from the 
zero-point vibration. This is given as 



A£(0*) = 



y/3a^N 



dkxdky 



IjItO- 7lstBZ 



(30) 



where a is the NN bond length, and J^^^g^.'^^xdky represents 
the integral in the 1st Brillouin zone. Numerical results of 
AE{9*) as a function of 9* are presented in Fig. [15] The en- 
ergy correction takes its minimum at six angles of 9* = nn/3 
with integer number «, reflecting the Ce symmetry in the 
free energyi^ It is worth noting that these are the same an- 
gle where the classical free energy takes the minimum (see 
Fig. |9l). That is, both the quantum and thermal fluctuations 
stabilize the same orbital configurations within the uniform PS 
alignments. Stability at these angles in the quantum model is 

attributed to the dispersion relation of the orbitons ©^ {9*); 
at 9* = nn/3, there is a one-dimensional zero-energy mode. 
For example, (0^-\9* = 0) = along the {k„ky) = (0,0) to 
(1,0) direction. This low-lying excitation contributes to the 
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FIG. 16: (a) Ground state energy and (b) energy gap for several 
size clusters. For the 2x3x3 cluster, energy difference between 
the doubly degenerate ground states and the first excited states are 
plotted. Broken lines are obtained by the least square fittings. 



energy gain from the quantum zero-point fluctuation. We sup- 
pose that, when the higher-order terms corresponding to the 
orbiton-obtion interactions, are taken into account, the disper- 
sion becomes gap-full, and the energy gain due to the zero- 
point fluctuation is reduced. 



B. Lanczos method 

In the numerical calculation in a finite size, the orbital struc- 
ture can be examined without assumption of the uniform PS 
configurations. In the quantum Monte-Cako simulation, we 
met a serious negative sign problem. Here we use the exact di- 
agonalization technique based on the Lanczos algorithm. We 
adopt finite size clusters of A^ = 2x2x2, 2x2x3, 2x3x3, 
and 2x3x4 sites with the periodic boundary condition. Be- 
cause of no conserved quantities in the Hamiltonian, all state 
vectors in the Hilbert space of 2^ dimension are dealt with in 
the Lanczos calculation. 

First, we show the ground state energy Eqs and the en- 
ergy gap A for several size clusters in Fig. [16] The ground 
state energy tends to approach, in the thermodynamic limit, 
around — 0.215A^y which is higher a little than the spin-wave 
results Eo + AE{e*) = -0.225NJ at 0* = n7t/3. Except for 
the 2x3x3 cluster, the ground state is not degenerate. The 
gap energy is defined as an energy difference between the 
ground state and the 1st excited one. The numerical value 
monotonically decreases with the system size A^, and seems 
to vanish in the thermodynamic limit. However, we cannot 
distinguish the two possibilities in an infinite system: degen- 
erate ground states and a non-degenerate one with gap-less 
excitation. The correlation functions of PS defined in Eq. ( |20] | 



1.0 



0.5 



0.0 
1.0 



0.5 



0.0 
1.0 



'0.5 



0.0 



(a) 



2x2x2 

(0,0) 



■ S™(10 
-1 S"(k) 



,0 



1 1 

2' 2 



^A 



(b) 2x3x3 



-5-5) (i -5) (o,o)(i°) {"-l) 






Hi) m 



(c) 2x3x4 



-.0 ±.0 0,- 



3 2I Vs 2 






FIG. 17: Correlation functions 5'"'(k) for several momenta k. Clus- 
ter sizes are 2 X 2 X 2 in (a), 2 x 3 x 3 in (b), and 2 x 3 x 4 in (c). 





FIG. 18: Some of the PS configurations where the honeycomb lattice 
is covered by NN bonds with the minimum bond energy. One of 
the ^ = 1 states in (a), and one of the q = —I in (b). In NN bonds 
surrounded by ellipses, the bond energy is the lowest. 



are calculated for several momenta and system sizes (Fig.fTTl). 
In the smallest size of 2 x 2 x 2 sites, 5""(-"' (k) at k = (0,0) 
stands out. However, with increasing A^, 5''"'(k)'s become al- 
most momentum independent and all of the values are less 
than 25% of the maximum. Reduction of S''" (k) at k = (0, 0) 
is faster than 1 /N. Thus, the conventional long-range order 
characterized by the correlation functions does not exist, as 
we have shown in the classical model. 

In the quantum system, the operator corresponding to q = 
^' Li cos 30, defined in Eq. (fTST i becomes a constant C- 
number due to the algebra for the 5=1/2 spin operator. Then, 
we adopt the variational-like method to analyze the ground- 
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FIG. 19: One example for the two PS configurations where a reso- 
nance state is possible due to the off-diagonal matrix element. 
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State wave function. As explained in Sect. Hill the classical PS 
states below Tq is characterized by the parameter q defined in 
Eq. (ISTT i. i.e. the PS angles are fixed at 2«7r/3 or (2« + l)7r/3 
with an integer number n. From these results and snapshots 
of the MC simulation, we consider the trial PS states where a 
honeycomb lattice is covered by the NN bonds with the mini- 
mum bond energy. Some examples are shown in Fig. [18] We 
construct the wave function as a linear combination of these 
states. This is given by 

|>p(±))=^£^,||V/(T))±|V/(^))|, (31) 

where -jV is a normalized factor, ^ are variational parame- 

ters, and | Vv ) is the wave function for the /-th PS configu- 
ration which satisfy the condition explained above. The wave 

function 1 1//| ') is given by the unitary transformation from the 
all-up PS state | f • • • f) as follows. 



IW 



,a)\ 



(T)\ 



n^('^'j)o7)jT •••!). 



(32) 



{ij)l 



Similarly, 1 1//,'^^') is obtained from the all-down state | i • • • i 
). The a bond direction is taken as the quantized axis, and 
the subscript {ij)[ represents the NN ij pair in the Z-the PS 
configuration. The operator U{^r\)iij), describes a rotation of 
T, and T; with respect to the T^' axis defined by 



U{'^n){ij),=^^V 



Tl 



Tj 



(33) 



where rj indicates a direction connecting / and j, and 
{^aA^^^y) == (0,2;r/3,47r/3). Because of the off-diagonal 

matrix elements among some states in |vv ) ^^^ !•// )' ^^^' 
tain kinds of resonance states are realized. A set of two PS 
configurations, termed 1^//,) and | <//«), shown in Fig. [T9l is 
an example. The off-diagonal matrix element between the 
two is (v/L|^|rff) == -yA^3/(16-2*'). This is about 10% 
of the energy gain due to quantum effect, (£gs — Eq)/JN, 
where /sgs is the ground state energy shown in Fig. [T6l a) and 
£0 = —3JN/16. Since the Hamiltonian is invariant under the 
inversion of all PS operators, the energy eigen states are classi- 
fied by the parity of this operation. The wave functions |^*+') 
and 1^'^') have the even and odd parities, respectively. Ex- 
cept for the degenerate ground state in the 2x3x3 size clus- 
ter, the ground state wave function |0) shows the even parity. 



FIG. 20: Overlap integrals between the ground-state wave function 
and the trial functions. Squares are for the results where the varia- 
tional parameters ia?/ are taken to be one, and circles are obtained by 
optimizing iS^. 



The doubly degenerate ground states in 2 x 3 x 3 are classi- 
fied as the even and odd parity states, and the even-parity one 
is used for the analyses. Figure |20] shows the overlap inte- 
gral IV ee |(0|1'W)|2 as a function of 1/Af. in a 2 x 2 x 2 size 
cluster, the ground-state wave function is almost completely 
reproduced by the trial function. With increasing A^, a value 
of W is gradually reduced. However, this reduction is rather 
weak by optimizing the variational parameters £/[, and W is 
maintained around 0.8 even in the largest size cluster. Thus, 
at least within the present calculation, the ground-state wave 
function is well reproduced by the trial wave function where 
the honeycomb lattice is covered by NN bonds with the mini- 
mum bond energy. 



V. DISCUSSION AND SUMMARY 

First, we have some remarks on low temperature orbital 
state in the classical model. As shown in Sect. [Till below 
To, the PS angle at each site is fixed at one of the three an- 
gles 2«;r/3 or one of the three (2« + l)7r/3. Here we discuss 
whether all (7 = ± 1 states appear equivalently or some specific 
PS configurations in the ^ = ± 1 states are stabilized. First we 
are able to exclude a possibility of the so-called directional 
order (DO). This is well examined in the orbital compass 



■ 30.31.32 



one com- 



model in a two-dimensional square lattice 
ponent in the PS operator, e.g. T', is aligned uniformly in 
each one-dimensional chain along a direction in the square 
lattice, e.g. the z direction, but there is no PS correlation be- 
tween the different chains. A natural order parameter of DO 



is a 



compass 



Li{TiT^,, 



T,^T-l^] ■ Below the DO temper- 
ature, a PS-angle parameter^, cos 20,, such as q in Eq. ( 1211 1. is 
developed, but the conventional PS correlation functions, such 
as ^'"'(k) in Eq. ( l20l i. are not. We introduce the honeycomb 
lattice version of the directional order parameter: 



D 



ieA 17 



t,''t/'. 



J2n^n/3 



(34) 
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When PS's are aligned uniformly inside zigzag chains along 
rj direction, but these chains are independent with each other, 
D acts as a monitor However, calculated (D^)''^ by the MC 
method are less than 5 x 10^^ and quickly disappears with 
increasing the system size. We also consider that a possibility 
of the Kosterlitz-Thouless transition^^*^ at To is weak. We 
calculate the uniform susceptibility 



Xu 



TN 



L(TrT, 



(35) 



and the corresponding correlation length ^„. Though the 
present system size is limited up to 2 x 24 x 24 sites, both 
Xu and ^„ do not show anomalous behavior around To, and 
their values decrease with increasing the system size. 

Here, we suggest that some topological PS configurations 
in hexagons may be more stabilized than other ^ = ± 1 states 
in the classical model. In a MC snapshot (see Fig. 12111. we of- 
ten find two characteristic PS configuration patterns in a hon- 
eycomb lattice; a uniform PS array termed configuration I, 
and a regular array of hexagons with maximum energy gain, 
termed configuration 11. Of course, their simple long-range 
orders are excluded from the calculated results of S''" (k) in 
Fig-El There are possibilities that two configurations coexist, 
and/or these are distributed randomly. These are monitored by 
a parameter «min which represents a number of NN bonds with 
the minimum bond energy in a hexagon. In the configuration 
II, hexagons with rimin = 3 and are aligned regularly. It is 
convenient to introduce the following parameter defined in a 
hexagon at r; 



Rir) 



9 (2 



8 3 



'-N{r) - 1 



1 



with 



A^(r)=I 



('■;■) 



3 '■ J 3 



(36) 



(37) 



where Y.(ij) represents a sum for six NN bonds in a hexagon. 
Because (t^- '^ T ■ ''' ) = 1/4 when a NN ij bond takes the min- 
imum bond energy, we have {N{r)) = «min- The parameter 
R{r) takes one for the hexagons with n^i„ = 0, 3, and zero for 
the hexagon with «rnin = 1,2. We calculate A^^'(Lr^(r)) in 
the 2x9x9 site cluster by the MC method. The calculated 
value is about 0.42 below Tq which is larger than a value (0.3) 
in the states where all ^ = ± 1 configurations appear equiva- 
lently. That is, the configuration II is expected to be more 
stabilized than other ^ = ±1 states. This is due to their low- 
energy fluctuations. Hexagons characterized as «niin = are 
included in the configuration II. As shown in Fig.|22l there are 
two-ways of fluctuation in each hexagon with «„,!„ = 0. When 
we consider the configuration II containing m hexagons with 
«min = 0, a number of configurations are roughly ;v/6Cm2'". 
This is remarkable in comparison with that in the configura- 
tion I; as explained in Fig. [8] there are also two ways of fluctu- 
ation in each zigzag chain. This corresponds to the so-called 
stacking degeneracy observed in the Cg orbital modeli^ When 
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FIG. 21: A MC snapshot for PS configurations. Numbers in each 
hexagons denotes a number of the minimum energy bonds «,nin- 




FIG. 22: Zero-dimensional fluctuation in the configuration II. Left: 
PS configuration. Right: a configuration obtained by ±59 rotation 
of PS's in each hexagon with «,„;„ = 0. 



we consider the configuration I containing m zigzag chains, 
a number of configurations are roughly /^C,„2"'. Difference 
between the two is attributed to dimensionality of the fluctu- 
ations. This zero-dimensional fluctuation seen in the config- 
uration II is unique in this honeycomb lattice model, and is 
expected to be an origin of no conventional long-range order. 
In summary, we study the doubly degenerate orbital model 
on a honeycomb lattice, motivated from an orbital state in 
multiferroic layered iron oxides /?Fe204. There is a macro- 
scopic number of degeneracy in the classical ground state, as 
seen in the three-dimensional eg orbital model. We mainly 
focus on lifting of the degeneracy due to thermal and quan- 
tum effects. In the classical and quantum spin-wave anal- 
yses, where the uniform orbital configurations are assumed, 
results are similar to those in the e^j-orbital model. Both ther- 
mal and quantum fluctuations stabilize the states with the PS 
angles of 9* = nn/3. Beyond the uniform configuration as- 
sumption, we apply the Monte-Carlo simulation to the clas- 
sical model. A peak structure in the specific heat is found 
around To /J = 0.006. However, below 7b, the PS correla- 
tion functions indexed by any possible momenta in clusters 
are not developed, unlike those in the e,, -orbital model. We 
find that the correlation function of a parameter for the orbital 
PS angle, q = ^, cos 30,, grows up below 7b, and reaches its 
maximum at the low temperature limit. That is, the PS an- 
gle at each site takes one of the three angles 2nn/3 or one of 
the three 2(« + l);r/3. This degeneracy lifting is attributed 
to existence of low-lying fluctuation around these configu- 
rations. We suggest that zero-dimensional fluctuation in the 



12 



hexagons plays a crucial role to make difference between the 
present honeycomb-lattice model and the eg orbital one in a 
cubic lattice. We also analyze the quantum model by utiliz- 
ing the Lanczos method. As seen in the classical model, any 
remarkable features are not shown in the two-body PS cor- 
relation functions. This suggests no conventional long-range 
order indexed by specific momenta. The ground-state wave 
function is analyzed by a variational-like method. This is well 
represented by a linear combination of the wave functions 
where a honeycomb lattice is covered by NN dimers with the 
minimum-energy PS configurations. 
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